F distribution (f)#
The F (Fisher–Snedecor) distribution in SciPy (scipy.stats.f) is a continuous distribution on \((0,\infty)\) that arises as a ratio of two independent scaled chi-square variables.
It is the workhorse behind ANOVA, overall regression significance tests, and classic tests comparing variances.
Learning goals#
Understand the chi-square ratio generative story and why the distribution is heavy-tailed.
Write down the PDF/CDF (and the Beta-function connection).
Derive the mean/variance and understand when moments do not exist.
Sample from \(\mathrm{F}(\mathrm{dfn},\mathrm{dfd})\) using a NumPy-only algorithm.
Use
scipy.stats.ffor evaluation, simulation, and fitting.
import numpy as np
import plotly
import plotly.express as px
import plotly.graph_objects as go
import os
import plotly.io as pio
import scipy
from scipy import optimize
from scipy import stats
from scipy.special import beta as beta_func
from scipy.special import betainc, betaln, digamma
pio.templates.default = "plotly_white"
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")
np.set_printoptions(precision=5, suppress=True)
rng = np.random.default_rng(42)
# Record versions for reproducibility (useful when numerical details matter).
VERSIONS = {"numpy": np.__version__, "scipy": scipy.__version__, "plotly": plotly.__version__}
VERSIONS
{'numpy': '1.26.2', 'scipy': '1.15.0', 'plotly': '6.5.2'}
1) Title & Classification#
Name:
f(F distribution; Fisher–Snedecor; SciPy:scipy.stats.f)Type: Continuous
Support (standard form): \(x \in (0,\infty)\)
Parameter space (standard form): numerator df \(\mathrm{dfn} > 0\), denominator df \(\mathrm{dfd} > 0\)
SciPy location/scale:
loc \in \mathbb{R},scale > 0with $\(X = \text{loc} + \text{scale}\,Y, \qquad Y \sim \mathrm{F}(\mathrm{dfn},\mathrm{dfd}).\)$
Unless stated otherwise, this notebook works with the standard form (loc=0, scale=1).
2) Intuition & Motivation#
What it models#
The F distribution models a ratio of (noise) energy levels.
A standard generative story is:
Draw \(U \sim \chi^2_{\mathrm{dfn}}\) and \(V \sim \chi^2_{\mathrm{dfd}}\) independently.
Normalize each by its degrees of freedom.
Take their ratio:
If \(U/\mathrm{dfn}\) and \(V/\mathrm{dfd}\) are thought of as variance-like quantities (they are averages of squared standard normals), then \(X\) is a variance ratio.
Typical real-world use cases#
ANOVA: compares explained variance to unexplained variance.
Regression F-test: tests whether a set of predictors improves fit beyond a baseline model.
Comparing variances: under Gaussian assumptions, \(s_1^2/s_2^2\) follows an F law under the null \(\sigma_1^2=\sigma_2^2\).
Relations to other distributions#
t distribution: if \(T \sim t_{\nu}\) then \(T^2 \sim \mathrm{F}(1,\nu)\).
Reciprocal symmetry: if \(X \sim \mathrm{F}(d_1,d_2)\) then \(1/X \sim \mathrm{F}(d_2,d_1)\).
Beta prime / Beta connection: if \(B \sim \mathrm{Beta}(a,b)\) with \(a=d_1/2\), \(b=d_2/2\), then $\(X = \frac{d_2}{d_1}\,\frac{B}{1-B} \sim \mathrm{F}(d_1,d_2).\)$ This identity explains why the CDF involves the regularized incomplete beta function.
3) Formal Definition#
Let \(d_1=\mathrm{dfn}>0\) and \(d_2=\mathrm{dfd}>0\).
PDF#
For \(x>0\) the F distribution has density
where \(\mathrm{B}(\cdot,\cdot)\) is the Beta function.
CDF#
Define the transformation
Then the CDF is
where \(\mathrm{I}_z(a,b)\) is the regularized incomplete beta function.
def f_logpdf(x: np.ndarray, dfn: float, dfd: float) -> np.ndarray:
"""Log-PDF of the standard F(dfn, dfd) distribution (loc=0, scale=1)."""
dfn = float(dfn)
dfd = float(dfd)
if dfn <= 0 or dfd <= 0:
raise ValueError("dfn and dfd must be > 0")
x = np.asarray(x, dtype=float)
out = np.full_like(x, -np.inf, dtype=float)
mask = x > 0
xm = x[mask]
a = 0.5 * dfn
b = 0.5 * dfd
ratio = dfn / dfd
# log f(x) = a log(dfn/dfd) - log B(a,b) + (a-1) log x - (a+b) log(1 + (dfn/dfd) x)
out[mask] = (
a * np.log(ratio)
- betaln(a, b)
+ (a - 1.0) * np.log(xm)
- (a + b) * np.log1p(ratio * xm)
)
return out
def f_pdf(x: np.ndarray, dfn: float, dfd: float) -> np.ndarray:
"""PDF of the standard F(dfn, dfd) distribution (loc=0, scale=1)."""
return np.exp(f_logpdf(x, dfn, dfd))
def f_cdf(x: np.ndarray, dfn: float, dfd: float) -> np.ndarray:
"""CDF of the standard F(dfn, dfd) distribution (loc=0, scale=1)."""
dfn = float(dfn)
dfd = float(dfd)
if dfn <= 0 or dfd <= 0:
raise ValueError("dfn and dfd must be > 0")
x = np.asarray(x, dtype=float)
out = np.zeros_like(x, dtype=float)
mask = x > 0
xm = x[mask]
a = 0.5 * dfn
b = 0.5 * dfd
# Stable transform u = dfn*x/(dfn*x+dfd) = 1 / (1 + dfd/(dfn*x))
u = 1.0 / (1.0 + (dfd / (dfn * xm)))
out[mask] = betainc(a, b, u)
return out
# Sanity check: our formulas match SciPy.
dfn, dfd = 5.0, 10.0
x = np.logspace(-3, 2, 25)
dist = stats.f(dfn, dfd)
assert np.allclose(f_pdf(x, dfn, dfd), dist.pdf(x))
assert np.allclose(f_cdf(x, dfn, dfd), dist.cdf(x))
assert np.allclose(f_logpdf(x, dfn, dfd), dist.logpdf(x))
4) Moments & Properties#
A key fact is that the F distribution has polynomial tails, so not all moments exist.
Let \(a=d_1/2\) and \(b=d_2/2\).
Raw moments#
For \(k < b\) (equivalently \(d_2 > 2k\)), the \(k\)-th raw moment exists and is
This immediately implies:
mean exists iff \(d_2>2\)
variance exists iff \(d_2>4\)
skewness exists iff \(d_2>6\)
(excess) kurtosis exists iff \(d_2>8\)
Mean / variance / skewness / kurtosis#
When the corresponding moments exist:
where \(\gamma_2\) is the excess kurtosis.
Mode#
If \(d_1>2\), the mode is
MGF / characteristic function#
The moment generating function \(M_X(t)=\mathbb{E}[e^{tX}]\) does not exist for any \(t>0\) (the tail is too heavy).
The Laplace transform \(\mathbb{E}[e^{tX}]\) exists for \(t<0\).
The characteristic function \(\varphi_X(t)=\mathbb{E}[e^{itX}]\) exists for all real \(t\) because \(|e^{itX}|\le 1\).
A generic integral representation is
Closed forms involve special functions (hypergeometric functions) and are rarely used directly in practice.
Entropy#
The differential entropy has a closed form:
where \(\psi\) is the digamma function.
def f_moment_raw(dfn: float, dfd: float, k: float) -> float:
"""Raw moment E[X^k] for X ~ F(dfn, dfd), when it exists (dfd > 2k)."""
dfn = float(dfn)
dfd = float(dfd)
k = float(k)
if dfn <= 0 or dfd <= 0:
raise ValueError("dfn and dfd must be > 0")
if dfd <= 2 * k:
return np.nan
a = 0.5 * dfn
b = 0.5 * dfd
return (dfd / dfn) ** k * np.exp(betaln(a + k, b - k) - betaln(a, b))
def f_mean(dfn: float, dfd: float) -> float:
return np.nan if dfd <= 2 else dfd / (dfd - 2)
def f_var(dfn: float, dfd: float) -> float:
if dfd <= 4:
return np.nan
return (2 * dfd**2 * (dfn + dfd - 2)) / (dfn * (dfd - 2) ** 2 * (dfd - 4))
def f_skew(dfn: float, dfd: float) -> float:
if dfd <= 6:
return np.nan
return ((2 * dfn + dfd - 2) * np.sqrt(8 * (dfd - 4))) / ((dfd - 6) * np.sqrt(dfn * (dfn + dfd - 2)))
def f_kurt_excess(dfn: float, dfd: float) -> float:
if dfd <= 8:
return np.nan
num = 12 * (dfn * (5 * dfd - 22) * (dfn + dfd - 2) + (dfd - 4) * (dfd - 2) ** 2)
den = dfn * (dfd - 6) * (dfd - 8) * (dfn + dfd - 2)
return num / den
def f_entropy(dfn: float, dfd: float) -> float:
a = 0.5 * float(dfn)
b = 0.5 * float(dfd)
return (
np.log(beta_func(a, b))
- (a - 1) * digamma(a)
- (b + 1) * digamma(b)
+ (a + b) * digamma(a + b)
+ np.log(dfd / dfn)
)
dfn, dfd = 7.0, 25.0
scipy_m, scipy_v, scipy_s, scipy_k = stats.f(dfn, dfd).stats(moments="mvsk")
out = {
"mean": (f_mean(dfn, dfd), scipy_m),
"var": (f_var(dfn, dfd), scipy_v),
"skew": (f_skew(dfn, dfd), scipy_s),
"kurt_excess": (f_kurt_excess(dfn, dfd), scipy_k),
"entropy": (f_entropy(dfn, dfd), stats.f(dfn, dfd).entropy()),
}
out
{'mean': (1.0869565217391304, 1.0869565217391304),
'var': (0.4822344816943791, 0.4822344816943791),
'skew': (1.741779266684047, 1.741779266684047),
'kurt_excess': (5.791950464396285, 5.7919504643962885),
'entropy': (0.8571957656704647, 0.8571957656704612)}
5) Parameter Interpretation#
The parameters are degrees of freedom:
\(d_1=\mathrm{dfn}\) controls the numerator chi-square \(U \sim \chi^2_{d_1}\).
\(d_2=\mathrm{dfd}\) controls the denominator chi-square \(V \sim \chi^2_{d_2}\).
How shape changes#
Small \(d_2\) means a noisy denominator \(V/d_2\), which creates a very heavy right tail.
This is why the mean/variance/skewness/kurtosis stop existing as \(d_2\) crosses \(2,4,6,8\).
Larger \(d_1\) pushes the distribution away from 0 (the density behaves like \(x^{d_1/2-1}\) near 0).
Useful limits#
As \(d_2 \to \infty\), the denominator concentrates: \(V/d_2 \to 1\), so
\[\mathrm{F}(d_1,d_2) \Rightarrow \chi^2_{d_1}/d_1.\]As both \(d_1,d_2\) grow, both \(U/d_1\) and \(V/d_2\) concentrate near 1, and the ratio concentrates near 1.
Reciprocal symmetry: \(X \sim \mathrm{F}(d_1,d_2)\) implies \(1/X \sim \mathrm{F}(d_2,d_1)\).
x = np.linspace(0.001, 6.0, 600)
# Vary dfd (tail heaviness)
dfn_fixed = 5
params_tail = [(dfn_fixed, 5), (dfn_fixed, 10), (dfn_fixed, 30)]
fig = go.Figure()
for dfn, dfd in params_tail:
fig.add_trace(go.Scatter(x=x, y=f_pdf(x, dfn, dfd), mode="lines", name=f"dfn={dfn}, dfd={dfd}"))
fig.update_layout(
title="F PDF: effect of the denominator df (tail heaviness)",
xaxis_title="x",
yaxis_title="pdf(x)",
)
fig.show()
# Vary dfn (left/peak behavior)
dfd_fixed = 10
params_peak = [(2.5, dfd_fixed), (5, dfd_fixed), (20, dfd_fixed)]
fig = go.Figure()
for dfn, dfd in params_peak:
fig.add_trace(go.Scatter(x=x, y=f_pdf(x, dfn, dfd), mode="lines", name=f"dfn={dfn}, dfd={dfd}"))
fig.update_layout(
title="F PDF: effect of the numerator df (behavior near 0 and peak)",
xaxis_title="x",
yaxis_title="pdf(x)",
)
fig.show()
6) Derivations#
This section sketches three useful derivations: the CDF, the mean/variance, and the likelihood.
Deriving the CDF via a Beta transformation#
Let \(X \sim \mathrm{F}(d_1,d_2)\) and define
Using the Beta-prime identity (or a Jacobian change of variables), one can show
Therefore
Deriving moments#
From the Beta connection
Then for \(k<b\):
Setting \(k=1\) and \(k=2\) yields \(\mathbb{E}[X]\) and \(\mathbb{E}[X^2]\); then
Likelihood (i.i.d. sample)#
Given data \(x_1,\dots,x_n\) i.i.d. from \(\mathrm{F}(d_1,d_2)\) (standard form), the log-likelihood is
There is no closed-form MLE for \((d_1,d_2)\); it is typically found by numerical optimization.
def f_loglik(dfn: float, dfd: float, data: np.ndarray) -> float:
data = np.asarray(data, dtype=float)
if np.any(data <= 0):
return -np.inf
return float(np.sum(f_logpdf(data, dfn, dfd)))
# Demonstration: recover parameters from synthetic data using log-likelihood.
rng_fit = np.random.default_rng(123)
dfn_true, dfd_true = 6.0, 18.0
sample = stats.f(dfn_true, dfd_true).rvs(size=2000, random_state=rng_fit)
def neg_loglik(theta: np.ndarray) -> float:
# Optimize in log-space to enforce positivity.
dfn = float(np.exp(theta[0]))
dfd = float(np.exp(theta[1]))
return -f_loglik(dfn, dfd, sample)
theta0 = np.log([5.0, 10.0])
res = optimize.minimize(neg_loglik, theta0, method="Nelder-Mead")
dfn_hat, dfd_hat = np.exp(res.x)
# Compare to SciPy's built-in MLE (fix loc=0, scale=1).
dfn_fit, dfd_fit, loc_fit, scale_fit = stats.f.fit(sample, floc=0, fscale=1)
{
"true": (dfn_true, dfd_true),
"mle_manual": (dfn_hat, dfd_hat),
"mle_scipy_fit": (dfn_fit, dfd_fit),
"success": res.success,
}
{'true': (6.0, 18.0),
'mle_manual': (6.26112794339186, 16.819704450943895),
'mle_scipy_fit': (6.261003935389471, 16.820016818351448),
'success': True}
7) Sampling & Simulation#
The chi-square ratio representation gives a simple sampling recipe:
Sample \(U \sim \chi^2_{d_1}\) and \(V \sim \chi^2_{d_2}\) independently.
Return \(X = (U/d_1)/(V/d_2)\).
A chi-square random variable is a Gamma:
So sampling from F reduces to sampling from a Gamma distribution. Below we implement a NumPy-only Gamma sampler using the Marsaglia–Tsang method.
Marsaglia–Tsang (high level)#
For \(k\ge 1\) (shape):
Draw \(Z\sim\mathcal{N}(0,1)\) and set \(V=(1+cZ)^3\).
Accept/reject using a cheap test first, then a log test.
For \(k<1\), boost to \(k+1\) and apply a power transform.
This is a standard, fast algorithm used in many libraries.
def gamma_rvs_numpy(shape: float, size: int, rng: np.random.Generator, scale: float = 1.0) -> np.ndarray:
"""Sample Gamma(shape, scale) using NumPy only (Marsaglia–Tsang).
Parameters
----------
shape:
k > 0
size:
number of samples
rng:
NumPy Generator
scale:
theta > 0 (multiplicative scale)
"""
k = float(shape)
theta = float(scale)
if k <= 0:
raise ValueError("shape must be > 0")
if theta <= 0:
raise ValueError("scale must be > 0")
# k < 1: boost to k+1 and apply power transform
if k < 1:
g = gamma_rvs_numpy(k + 1.0, size, rng, scale=1.0)
u = rng.random(size)
return theta * g * (u ** (1.0 / k))
# k >= 1: Marsaglia–Tsang
d = k - 1.0 / 3.0
c = 1.0 / np.sqrt(9.0 * d)
out = np.empty(size, dtype=float)
filled = 0
while filled < size:
n = size - filled
x = rng.standard_normal(n)
v = 1.0 + c * x
v = v * v * v # (1 + c x)^3
u = rng.random(n)
positive = v > 0
# First (cheap) acceptance
accept = positive & (u < 1.0 - 0.0331 * (x**4))
# Second acceptance (log test)
log_v = np.zeros_like(v)
log_v[positive] = np.log(v[positive])
accept2 = positive & (~accept) & (np.log(u) < 0.5 * x * x + d * (1.0 - v + log_v))
accept = accept | accept2
accepted = d * v[accept]
take = min(accepted.size, n)
out[filled : filled + take] = accepted[:take]
filled += take
return theta * out
def chisquare_rvs_numpy(df: float, size: int, rng: np.random.Generator) -> np.ndarray:
"""Sample ChiSquare(df) using Gamma(df/2, scale=2)."""
df = float(df)
if df <= 0:
raise ValueError("df must be > 0")
return gamma_rvs_numpy(df / 2.0, size, rng, scale=2.0)
def f_rvs_numpy(dfn: float, dfd: float, size: int, rng: np.random.Generator) -> np.ndarray:
"""Sample F(dfn, dfd) via ratio of independent chi-squares."""
u = chisquare_rvs_numpy(dfn, size, rng)
v = chisquare_rvs_numpy(dfd, size, rng)
return (u / dfn) / (v / dfd)
# Monte Carlo validation against theory
n = 200_000
params = (5.0, 10.0)
samples_numpy = f_rvs_numpy(*params, n, rng)
m_theory, v_theory = stats.f(*params).stats(moments="mv")
np.mean(samples_numpy), np.var(samples_numpy), m_theory, v_theory
(1.2529035200365988, 1.3745254004774181, 1.25, 1.3541666666666667)
8) Visualization#
We will visualize:
the PDF for several parameter settings
the CDF and an empirical CDF from Monte Carlo samples
a histogram of simulated samples overlaid with the theoretical PDF
dfn, dfd = 5.0, 10.0
x = np.linspace(0.001, 6.0, 600)
# PDF line + Monte Carlo histogram
n_vis = 60_000
samples = f_rvs_numpy(dfn, dfd, n_vis, rng)
fig = go.Figure()
fig.add_trace(
go.Histogram(
x=samples,
histnorm="probability density",
nbinsx=120,
name="Monte Carlo",
opacity=0.55,
)
)
fig.add_trace(go.Scatter(x=x, y=f_pdf(x, dfn, dfd), mode="lines", name="PDF (theory)", line=dict(width=3)))
fig.update_layout(
title=f"F(dfn={dfn}, dfd={dfd}): histogram vs PDF",
xaxis_title="x",
yaxis_title="density",
bargap=0.02,
)
fig.show()
# CDF vs empirical CDF
xs = np.sort(samples)
ecdf = np.arange(1, xs.size + 1) / xs.size
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=f_cdf(x, dfn, dfd), mode="lines", name="CDF (theory)", line=dict(width=3)))
fig.add_trace(go.Scatter(x=xs, y=ecdf, mode="lines", name="CDF (empirical)", line=dict(width=2, dash="dot")))
fig.update_layout(
title=f"F(dfn={dfn}, dfd={dfd}): empirical CDF vs theory",
xaxis_title="x",
yaxis_title="CDF",
)
fig.show()
9) SciPy Integration (scipy.stats.f)#
SciPy provides a robust implementation with location/scale support.
from scipy import stats
dist = stats.f(dfn, dfd) # standard form
pdf_vals = dist.pdf(x)
cdf_vals = dist.cdf(x)
rv = dist.rvs(size=1000, random_state=rng)
# Fit (MLE). By default SciPy fits dfn, dfd, loc, scale.
# If you want the *standard* F family, fix loc=0 and scale=1.
dfn_hat, dfd_hat, loc_hat, scale_hat = stats.f.fit(data, floc=0, fscale=1)
Below is a short, end-to-end example.
dfn, dfd = 4.0, 12.0
x0 = np.array([0.25, 0.5, 1.0, 2.0, 4.0])
dist = stats.f(dfn, dfd)
pdf = dist.pdf(x0)
cdf = dist.cdf(x0)
rv = dist.rvs(size=5, random_state=rng)
# Fitting: generate data, then recover parameters.
rng_fit = np.random.default_rng(7)
data = stats.f(dfn, dfd).rvs(size=4000, random_state=rng_fit)
dfn_hat, dfd_hat, loc_hat, scale_hat = stats.f.fit(data, floc=0, fscale=1)
{
"x": x0,
"pdf": pdf,
"cdf": cdf,
"rvs": rv,
"fit": {"dfn_hat": dfn_hat, "dfd_hat": dfd_hat, "loc": loc_hat, "scale": scale_hat},
}
{'x': array([0.25, 0.5 , 1. , 2. , 4. ]),
'pdf': array([0.61496, 0.67983, 0.46719, 0.15676, 0.02124]),
'cdf': array([0.09586, 0.26351, 0.55505, 0.84137, 0.97256]),
'rvs': array([1.41411, 1.15258, 0.53981, 1.69128, 0.04488]),
'fit': {'dfn_hat': 4.006050540012662,
'dfd_hat': 11.86095144322162,
'loc': 0,
'scale': 1}}
10) Statistical Use Cases#
A) Hypothesis testing#
Variance ratio test (Gaussian assumption)
If \(X_1,\dots,X_{n_1}\) and \(Y_1,\dots,Y_{n_2}\) are independent normal samples with equal variance under \(H_0\), then
where \(S^2\) denotes the sample variance.
ANOVA / regression
F statistics often take the form
and compare explained vs unexplained variability.
B) Bayesian modeling (variance ratios)#
For a normal sample with Jeffreys prior \(p(\mu,\sigma^2)\propto 1/\sigma^2\), the posterior for \(\sigma^2\) is an inverse-chi-square. For two independent groups, the posterior for the variance ratio satisfies
This yields a closed-form credible interval for \(\sigma_1^2/\sigma_2^2\).
C) Generative modeling#
The F distribution is a flexible heavy-tailed positive distribution. One trick is to use it as a random scale:
This produces occasional large scales (outliers) while remaining easy to simulate.
# A) Variance ratio test example
rng_test = np.random.default_rng(0)
n1, n2 = 25, 18
sigma = 2.0
x = rng_test.normal(0.0, sigma, size=n1)
y = rng_test.normal(0.0, sigma, size=n2)
s1 = np.var(x, ddof=1)
s2 = np.var(y, ddof=1)
F_stat = s1 / s2
# Two-sided p-value
p_left = stats.f.cdf(F_stat, n1 - 1, n2 - 1)
p_right = stats.f.sf(F_stat, n1 - 1, n2 - 1)
p_two_sided = 2 * min(p_left, p_right)
# B) Bayesian variance ratio posterior (Jeffreys prior)
nu1, nu2 = n1 - 1, n2 - 1
scale_ratio = s1 / s2
# Posterior: (sigma1^2/sigma2^2) | data ~ scale_ratio * F(nu2, nu1)
ci_level = 0.95
alpha = 0.5 * (1 - ci_level)
q_lo = stats.f.ppf(alpha, nu2, nu1)
q_hi = stats.f.ppf(1 - alpha, nu2, nu1)
ci_post = (scale_ratio * q_lo, scale_ratio * q_hi)
# Monte Carlo check: sample sigma_i^2 ~ Inv-chi-square(nu_i, s_i^2)
mc = 80_000
sigma1_sq = (nu1 * s1) / chisquare_rvs_numpy(nu1, mc, rng_test)
sigma2_sq = (nu2 * s2) / chisquare_rvs_numpy(nu2, mc, rng_test)
ratio_mc = sigma1_sq / sigma2_sq
ci_mc = (np.quantile(ratio_mc, alpha), np.quantile(ratio_mc, 1 - alpha))
# C) Generative modeling: heavy-tailed scale mixture
m = 120_000
s = f_rvs_numpy(5.0, 5.0, m, rng_test)
eps = rng_test.standard_normal(m)
y_heavy = np.sqrt(s) * eps
y_normal = eps
# Compare tail probabilities
thr = 3.0
p_tail_heavy = np.mean(np.abs(y_heavy) > thr)
p_tail_normal = np.mean(np.abs(y_normal) > thr)
{
"variance_ratio_test": {"F": F_stat, "p_two_sided": p_two_sided},
"posterior_CI_closed_form": ci_post,
"posterior_CI_monte_carlo": ci_mc,
"tail_P(|Y|>3)": {"heavy": p_tail_heavy, "normal": p_tail_normal},
}
{'variance_ratio_test': {'F': 1.06647032697333,
'p_two_sided': 0.9077897055189864},
'posterior_CI_closed_form': (0.4166185566400728, 2.545110186554938),
'posterior_CI_monte_carlo': (0.418158015120004, 2.5405681465872147),
'tail_P(|Y|>3)': {'heavy': 0.033525, 'normal': 0.0026333333333333334}}
11) Pitfalls#
Invalid parameters: \(d_1\le 0\) or \(d_2\le 0\) is invalid.
Non-existent moments: if \(d_2\le 2\), the mean is infinite; if \(d_2\le 4\), the variance is infinite; etc.
Numerical stability: for extreme \(x\) or degrees of freedom, prefer
logpdfoverpdf.loc/scaleconfusion in SciPy:scipy.stats.fis a four-parameter family by default; fixloc=0, scale=1if you mean the standard F distribution.Fitting: MLE can be sensitive to outliers (heavy tail) and may converge slowly; check diagnostics and consider robust alternatives when appropriate.
12) Summary#
The F distribution models a ratio of two variance-like quantities: \((\chi^2_{d_1}/d_1)/(\chi^2_{d_2}/d_2)\).
The PDF/CDF have closed forms involving the Beta function and the regularized incomplete beta.
Moments exist only up to order \(k<d_2/2\); in particular, the mean requires \(d_2>2\) and the variance requires \(d_2>4\).
Sampling is straightforward via Gamma/chi-square sampling; a NumPy-only Marsaglia–Tsang Gamma sampler works well.
In practice,
scipy.stats.fprovides evaluation, simulation, and MLE fitting.
References
Marsaglia, G. & Tsang, W. W. (2000). A Simple Method for Generating Gamma Variables.
Johnson, Kotz, and Balakrishnan. Continuous Univariate Distributions, Volume 2 (2nd ed.), Wiley, 1995.
SciPy documentation:
scipy.stats.f,scipy.special.betainc.